NASA  Contractor  Report  198205 
ICASE  Report  No.  95-59 


'X 


ICASE 


PSEUDO-TIME  METHOD  FOR  OPTIMAL  SHAPE 
DESIGN  USING  THE  EULER  EQUATIONS 


Angelo  lollo 
Geojoe  Kuruvila 
Shlomo  Ta’asan 

ilppioved  toT  puDi:.c  r'^jc.r;- ;  '• 


NASA  Contract  No.  NAS 1-19480 
August  1995 

Institute  for  Computer  Applications  in  Science  and  Engineering 
NASA  Langley  Research  Center 
Hampton,  VA  23681-0001 

Operated  by  Universities  Space  Research  Association 


19951013  017 

National  Aeronautics  and 
Space  Administration 


DTiarmAr^ 


Langley  Research  Center 

Hampton,  Virginia  23681-0001 


Pseudo-Time  Method  for  Optimal  Shape  Design 
Using  the  Euler  Equations 


Angelo  lollo^ 

Dipartimento  di  Ingegneria  Aeronautica  e  Spaziale 
Politecnico  di  Torino 
10129  Torino,  Italy 

Geo  joe  Kuruvila^ 

Advanced  Transport  Aircraft  Development 
McDonnell  Douglas  Aerospace 
Long  Beach,  CA  90810  1870 

Shlomo  Ta’asan^ 

Department  of  Mathematics 
Carnegie  Mellon  University 
Pittsburgh,  PA  15213 


Accesion  For 


NTiS  CRA&I 
OTIC  TAB 

n 

Unannounced 

Justification 

□ 

By 

Distribution  / 

Availability 

Codes 

Dist 


A'l 


Avail  and/or 
Special 


Abstract 

In  this  paper  we  exploit  a  novel  idea  for  the  optimization  of  flows  governed  by  the 
Euler  equations.  The  algorithm  consists  of  marching  on  the  design  hypersurface  while 
improving  the  distance  to  the  state  and  costate  hypersurfaces.  We  consider  the  problem 
of  matching  the  pressure  distribution  to  a  desired  one,  subject  to  the  Euler  equations, 
both  for  subsonic  and  supersonic  flows.  The  rate  of  convergence  to  the  minimum  for 
the  cases  considered  is  3  to  4  times  slower  than  that  of  the  analysis  problem.  Results 
are  given  for  Ringleb  flow  and  a  shockless  recompression  case. 
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1  Introduction 


In  recent  years  there  has  been  a  renewed  interest  in  optimal  design  in  fluid  mechanics. 
Faster  computers  and  reliable  numerical  simulations  make  feasible  some  of  the  aerodynamics 
optimal  design  problems  which  are  of  engineering  interest. 

The  statement  above  becomes  only  partially  true  when  we  consider  either  flows  governed 
by  Euler  or  Navier-Stokes  equations,  or  complicated  geometrical  configurations  with  many 
control  parameters.  For  such  cases,  shape  optimization  seems  to  be  still  not  practical  due 
to  extremely  time  consuming  computation.  The  present  work  aims  at  a  flexible  and  feasible 
approach  for  such  intensive  computational  problems  by  applying  a  novel  algorithm  proposed 
by  Ta’asan  in  [19]. 

The  problem  of  finding  a  shape  that  achieves  given  performance  has  been  attacked  by 
means  of  inverse  problem  formulations  [12], [2], [21], [5].  These  methods  have  in  common  the 
advantage  of  being  solved  at  the  same  cost  of  an  analysis  problem.  They  are  in  general  not 
extendable  to  three  dimensions.  Moreover  the  set  of  problems  that  can  be  solved  by  means 
of  inverse  design  is  limited. 

A  more  general  framework  is  to  consider  aerodynamics  design  problems  as  optimization 
problems.  From  the  mathematical  viewpoint  the  problem  is  to  find  U  such  that 

f  U^U 

\  S{U)  <  S{V)  WeU 

where  ZY  is  a  given  set  and  5  is  a  real-valued  functional  defined  on  U. 

Shape  design  optimization  problems  are  tightly  related  to  control  of  a  system  governed 
by  partial  differential  equations  where  the  controls  are  on  the  boundary.  Lions  set  in  [13] 
the  mathematical  framework  for  such  problems.  The  theory  is  concerned  mainly  with  linear 
systems  and  is  devoted  “(i)  to  obtain  necessary  (or  possibly  necessary  and  sufficient)  con¬ 
ditions  for  U  to  be  an  extremum  (or  minimum),  (ii)  to  study  the  structure  and  properties  of 
the  equations  expressing  these  conditions,  (iii)  to  obtain  constructive  algorithms  amenable 
to  numerical  computations  for  the  approximation  of  t/” . 

Pironneau  ([15], [16])  derived  an  adjoint  method  for  the  minimum-drag  problem  in  Stokes 
flows  and  subsequently  in  flows  governed  by  the  incompressible  Navier-Stokes  equations. 
Since  a  Navier-Stokes  solver  was  not  available,  some  solutions  were  obtained  using  simpler 
models;  see  Glowinski  and  Pironneau  [6]. 

For  the  Euler  equations  Jameson  proposed  in  [10]  an  adjoint  method  for  wing  design 
which  makes  use  of  conformal  mapping  to  control  the  shape  of  the  wing.  lollo,  Salas  and 
Ta’asan  [9]  studied  the  case  of  Euler  flows  with  embedded  shocks  for  a  one-dimensional 
case,  and  discussed  the  boundary  conditions  for  the  adjoint  equations.  At  the  shock  location 
it  was  shown  that  further  conditions  are  needed  for  the  adjoint  equation  to  be  well  posed. 
Subsequently,  lollo  and  Salas  extended  these  results  to  two-dimensional  flows,  and  presented 
computations  with  higher-order  spatial  accuracy  [8]. 

The  high  computational  cost  for  solving  optimization  problems  governed  by  fluid  dynamics 
equations  comes  from  several  sources.  The  first  is  the  cost  of  a  single  analysis,  which  for  the 
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Navier-Stokes  equations  in  three  dimensions  is  of  the  order  of  a  few  CRAY  hours.  Another 
source  is  the  fact  that  a  repeated  solution  of  the  flow  equations  may  be  required  for  gradient 
methods.  An  additional  significant  cost  may  arise  from  the  calculation  of  the  gradient  of 
the  functional. 

The  use  of  adjoint  methods  eliminates  the  unnecessary  cost  resulting  from  computation 
of  the  gradient,  and  is  much  more  efficient  compared  to  other  methods  including  finite- 
differences  and  sensitivity  analysis.  It  requires  the  computation  of  an  extra  system  of 
partial  differential  equations  (PDEs),  namely  the  costate  equations,  but  the  total  cost  for 
gradient  calculation  is  independent  of  the  number  of  design  variables.  A  comparison  study 
of  calculating  gradients  using  adjoint  methods  and  finite-difference  methods  was  done  by 
Beux  and  Dervieux  [3].  They  also  solved  pressure  reconstruction  problems  for  compressible 
internal  flows,  comparing  the  performances  of  several  algorithms.  Flows  with  embedded 
shocks  were  not  considered  in  this  work. 

The  adjoint  method,  being  an  efficient  method  for  calculating  the  gradient,  does  not  ad¬ 
dress  the  computational  expense  related  to  the  number  of  gradient  iterations  required  to 
reach  the  minimum.  In  general,  the  number  of  iterations  required  to  achieve  the  minimum 
grows  more  than  linearly  with  the  number  of  controls  used,  making  infeasible  design  prob¬ 
lems  in  three  dimensions  with  many  design  variables. 

Ta’asan  proposed  in  [18]  an  algorithm  to  reduce  the  cost  of  the  optimization  to  that  of  a 
single  analysis,  namely  the  one  shot  method.  The  idea  is  to  solve  the  flow  equations,  the 
costate  equations  and  the  optimality  condition  at  the  same  time.  The  main  idea  in  that 
algorithm  was  to  perform  the  optimization  iteration  on  coarse  grids  that  are  used  anyway 
in  the  multigrid  process.  Small  numbers  of  design  variables  were  considered  in  that  case. 

Ta’asan,  Kuruvila  and  Salas  [20]  applied  this  technique  to  a  potential  flow,  and  extended 
the  method  to  cases  of  moderate  numbers  of  design  variables.  Different  design  variables 
are  associated  with  different  grids  depending  on  the  smoothness  of  the  shape  functions 
associated  with  them,  and  are  updated  on  these  grids.  The  performance  of  this  algorithm 
was  practically  independent  of  the  number  of  design  variables. 

Arian  and  Ta’asan  [1]  extended  the  one  shot  method  to  infinite-dimensional  design  space. 
The  main  idea  of  the  method  was  to  construct  a  relaxation  that  smoothes  the  errors  in 
the  design  variables.  Application  to  control  problems  and  shape  design  problems  have 
demonstrated  solution  of  the  full  optimization  problem  in  a  cost  comparable  to  that  of  solving 
the  analysis  just  a  few  times,  independent  of  the  number  of  design  variables  (experiments 
using  up  to  128  design  variables  have  been  done). 

Beux  and  Dervieux  [4]  proposed  a  hierarchical  strategy  in  which  the  number  of  control 
parameters  is  progressively  increased  performing  a  multilevel  optimization  that  seems  to 
render  the  computational  cost  independent  from  the  number  of  control  parameters. 

The  drawbacks  of  the  one  shot  methods  are  their  programming  complexity  and  the  fact 
that  their  use  is  limited  to  multigrid  solvers.  This  was  the  motivation  for  the  study  of  a  new 
type  of  solution  strategy  for  optimization  problems  governed  by  PDEs  [19].  The  goal  was  to 
try  to  get  methods  that  solve  the  optimization  problem  in  a  cost  comparable  to  that  of  the 
analysis.  The  emphasize  was  on  simplicity  and  flexibility  to  work  in  existing  frameworks 


2 


which  do  not  necessarily  involve  multigrid  methods. 

The  main  observation  is  the  following.  The  solution  of  the  optimization  problem  lies  on 
the  intersection  of  the  state,  costate  and  design  hypersurfaces  (in  the  state,  costate  and 
design  spaces).  Gradient-based  methods  (including  adjoint  formulations)  can  be  viewed  as 
marching  along  the  intersection  of  the  state  and  costate  hypersurfaces.  This  is  an  expensive 
process  since  each  step  requires  the  solution  of  two  PDEs.  The  idea  of  the  pseudo-time 
method  was  to  perform  the  marching  on  the  design  hypersurface  while  improving  the  dis¬ 
tance  to  the  other  two.  The  cost  of  such  an  iteration  per  step  is  significantly  smaller  than 
that  of  gradient-based  methods.  Its  convergence  has  been  shown  by  Ta’asan  in  [19]  to  be 
independent  of  the  number  of  design  variables. 

In  the  present  paper  we  apply  the  pseudo-time  method  to  optimization  problem  using  the 
Euler  equations.  Using  this  method  the  cost  of  optimization  becomes  of  the  same  order  as 
that  of  analysis.  Moreover  the  algorithm  may  be  implemented  with  no  substantial  changes 
to  existing  codes.  Numerical  results  indicate  that  the  method  converges  at  a  rate  which  is 
independent  of  the  number  of  design  variables. 

2  Problem  statement 

The  Euler  equations  are  given  by 

U(  -|-  Fx  -|-  Gj,  =  0  (1) 

/  P  \  (  \  (  \ 

U  =  F  =  -P  +  G  = 

pv  puv  p  -|-  pv^ 

\  pe  j  V  «(pe  4-  p)  /  v{pe  -{■  p)  ) 

p  —  density 

u  =  x-component  of  the  velocity  vector 
V  =  y-component  of  the  velocity  vector 
e  =  specific  total  energy 
p  =  pressure 
a  —  speed  of  sound 
7  =  ratio  of  specific  heats 
7-  1 

^  =  — 

and  p  —  Kp{2e  —  v?  —  v^).  Furthermore  let 

au  =  au  = 


where 


with 
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(2) 


We  assume  that  these  equations  are  defined  on  a  domain  $  which  includes  a  sub-domain 
0,  whose  boundary  is  denoted  by  P.  On  the  boundary  we  define  a  curvilinear  coordinate  s 
and  a  normal  n  =  (nx,ny)  pointing  outward.  A  real  valued  functional  £(r, !/(?))  is  given, 
where  lf(r)  is  the  solution  of  the  Euler  equations  with  boundary  conditions  on  P.  The 
optimization  problem  that  we  study  is 

minimize  the  functional  £(P,y(P))  over  all  the  admissible  shapes  of  the  boundary  P. 

We  focus  on  the  following  model  problem.  The  sub-domain  0,  is  represented  by  a  nozzle; 
see  Fig.  1.  At  the  inlet,  total  pressure,  total  temperature  and  the  ratio  a  =  v/u  are  fixed. 
At  the  outlet,  if  the  flow  is  subsonic,  the  static  pressure  is  fixed  and  at  the  solid  walls  the 
impermeability  condition  unx  +  vny  =  0  is  enforced.  The  upper  wall  is  kept  fixed.  The 
lower  wall  0  is  represented  by  mean  of  the  parameterization 

J/(®)  =  (3) 

i 


where  the  functions  fi{x)  are  some  shape  functions  and  a  =  (ai, ...,  a,-, ...)  is  the  corres¬ 
ponding  set  of  shape  coefficients.  Given  a  desirable  lower  wall  pressure  distribution  p*{x) 
and  denoting  by  p®(x)  the  actual  one  on  the  lower  wall,  the  optimization  problem  consists 
in  finding  a  set  of  shape  coefficients  such  that  the  functional 

^  (4) 


is  minimized. 


3  Optimality  conditions 

The  optimality  conditions  are  derived  by  introducing  Lagrange  multipliers  and  considering 
the  augmented  functional 

£(U,  a,  A,  //)  =  S{a,  U)  +  /  '  A(  AU,  -b  BU^)dO  +  /  /xpV  •  nds  (5) 

Ja  JQ 

where  V  =  (u,  u).  The  vector  A{x,  y)  =  *(Ai,  A2,  A3,  A4),  and  the  scalar  fi{s)  are  the  Lagrange 
multipliers. 

Calculating  the  variation  of  the  functional  C  with  respect  to  the  variation  of  the  functions 
U,  A,  p,  and  the  parameters  a,-  respectively,  we  obtain  (see  [8]) 


5Cu 


(p®  —  p*)\5dx  -f-  f  ^A(Anx  -t-  Bny)Uds  -|- 

0 

/j'AxA  +  %B)UdQ  +  ^  An  ^  Uds 


A  dp 

Ja  W 


(6) 


4 


where 


^-2k( +  -u  -vl]  and  ^  °  ^  °  ^ 

dV~^^[  2  ’  5U  lO  0  1  0 


and 


5Ca  =  /  'A(AU^  +  BVy)dn 

J  ri 

<5£„  =  /  yupV  •  nds 

Je 

E  f  ^  (p®  -  P*)  fi  dx+  [  ‘A(AU,  +  BVy)  U  cos  6  ds-{- 

■  \Ja  dy  Q  Je 


5Ca  = 


+  /  P 


f  ■  nfids  —  f  ypY  ■  t  ^  cos^  6ds  + 

Je  oy  Je  dx 


dy 


+  /  ypY  -n^  sin^cos^ds  )  d,- 

Je  dx  ) 


(7) 

(8) 


(9) 


where  0  is  the  angle  between  the  normal  n  and  the  y-axis,  t  =  {—Uy,  rix),  and  U,  A,  y  and 
di  are  the  variations  of  the  corresponding  arguments.  _ 

At  the  minimum  of  the  functional,  for  all  the  possible  choices  of  the  functions  U,  A,  p 
and  of  the  parameters  d,  we  must  have 


5jC.u  =  SjC\  =  SHy  =  SjCoi  =  0. 


(10) 


Therefore,  we  have 
and 


6£a  =  0  ^  AUa;  +  BUy  =  0  on  f2 


5Cfi  =  0  4^  pY  •  n  =  0  on  0 

which  are  the  Euler  equations  and  boundary  conditions.  Furthermore 

5Cjj  =  0  ^ AAa;  +  ^BAy  =  0  on  Cl 

and 


dp 


where 


(p"'  -  p*)  cos  9  +  ^A(Ana:  +  Bny)  +  pn  =  0  on  0 


=  —  j^Ai  +  u\2  +  vXs  +  (qe  —  kV^)\4 


(11) 

(12) 

(13) 


For  the  boundary  condition  on  inlet  and  outlet  we  refer  the  reader  to  [8].  Given  U,  the  set 
of  costate  eqs. (11-13)  determine  uniquely  A  in  fl  and  y  on  0.  Finally,  given  a  and  knowing 
U  and  A,  we  can  calculate  from  eq.(9) 

dC  /■*’  dp 


dai 


[  ^  (P®  -  P*)  U  dx+  [  'A(AU,  +  BUy)  fi  cos  9  ds  + 

Ja  dy  0  7© 


(14) 


+  f  ■  nfids  —  f  ixpV  ■  t  ~  cos^  9ds  + 

JQ  oy  Je  dx 

+  /  upV  •  n  — ^  sin  0  cos  0  ds 

Je  dx 

In  case  of  a  shock  occuring  in  the  flow  field,  we  split  the  domain  of  integration  by  means  of 
a  curve  T  that  coincides  with  the  shock  where  it  exists.  Then  we  follow  the  same  derivation 
presented  so  far  on  each  of  the  two  sub-domains,  regarding  T  as  a  boundary;  see  [8].  The 
resulting  extra  condition  for  A  on  the  shock  is  A  =  0.  It  should  be  noted  that  if  the 
shocks  are  not  treated  properly,  the  problem  of  solving  the  costate  equations  with  boundary 
conditions  is  not  well-posed.  Jameson  presented  in  [11]  results  for  transonic  flows  over 
airfoils  where  the  wave-drag  is  minimized.  He  does  not  use  any  special  treatment  for  the 
shock  but  the  costate  equations  converge.  This  is  due  to  the  fact  that  the  scheme  that 
he  uses  for  solving  the  Euler  equations  smears  the  shocks  over  several  grid  points,  due  to 
artificial  viscosity. 


4  Pseudo- time  optimization  method 

There  are  many  methods  for  obtaining  the  minimum  of  the  functional  £  knowing  its  gradient 
with  respect  to  the  controls.  Adjoint  methods  involve  the  following  steps; 

1.  Start  with  a  set  a,-  of  shape  coefficients 

2.  Enforce  5£\  =  0  and  =  0  by  finding  a  U  that  satisfies  the  steady  state  Euler 
equations  and  boundary  conditions 

3.  Enforce  SCjj  =  0  by  finding  a  A  that  satisfies  the  costate  equations  and  boundary 
conditions 

4.  Calculate  Vo,£,  if  it  is  0  we  have  found  the  minimum,  otherwise 

5.  Update  a  with  Va£,  using  a  proper  stepsize. 

6.  Restart  from  2  until  Va£  =  0. 

The  need  to  repeat  steps  2  and  3  above  many  times  can  become  prohibitively  expensive 
for  geometrically  complex  configurations  requiring  computational  power  near  the  limits  of 
present  capabilities. 

Ta’cisan  proposed  in  [19]  an  efficient  way  of  solving  the  optimization  problem.  The  main 
observation  is  the  following.  The  solution  of  the  optimization  problem  lies  on  the  intersection 
of  the  state,  costate  and  design  hypersurfeces  (in  the  state,  costate  and  design  spaces). 
Gradient  based  methods  (including  adjoint  formulations)  can  be  viewed  as  marching  along 
the  intersection  of  the  state  and  costate  hypersurfaces.  This  is  an  expensive  process  since 
each  step  requires  the  solution  of  two  PDEs.  The  idea  of  the  pseudo-time  method  was  to 
perform  the  marching  on  the  design  space  while  improving  the  distance  to  the  other  two. 
Compare  Fig.  2  and  3.  The  cost  of  such  an  iteration  per  step  is  significantly  smaller  than 
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that  of  gradient  based  methods.  Its  convergence  has  been  shown  by  Ta’asan  in  [19]  to  be 
independent  of  the  number  of  design  variables. 

The  design  equation  for  a  wide  class  of  problems,  including  the  one  considered  here,  is 
defined  on  the  boundary  only.  Thus,  it  can  be  viewed  as  an  extra  boundary  condition, 
and  the  design  variables  as  the  additional  variables  to  solve  for.  In  some  cases  the  design 
equation  can  be  solved  for  the  design  variables  and  a  simple  implementation  of  the  above 
idea  exists.  In  other  cases  the  design  equation,  viewed  as  an  equation  for  the  design  variables 
keeping  the  state  and  costate  fixed,  may  be  singular  and  a  more  involved  implementation  is 
required.  This  is  the  case  for  the  problem  considered  here.  In  such  cases  it  is  necessary  to 
solve  for  the  design  variables  together  with  the  state  and  costate  variables  in  a  small  vicinity 
of  the  boundary  S. 

Thus,  at  each  step  of  computation  on  the  entire  field,  the  design  equation  is  satisfied 
together  with  the  boundary  conditions  for  the  state  and  costate  equations.  The  solution  on 
the  entire  field  (I  affects  the  result  of  the  optimization  on  S  through  the  values  of  U  and  A 
on  the  auxiliary  boundary  see  Fig.  4. 

The  algorithm  is  as  follows: 

1.  Start  with  a  tentative  set  of  o;,-. 

2.  March  in  time,  on  the  entire  field,  the  state  equation  a  few  steps. 

3.  March  in  time,  on  the  entire  field,  the  costate  equation  a  few  steps. 

4.  Solve  in  S  the  state  equation  with  its  boundary  conditions,  the  costate  equation  with 
its  boundary  conditions  and  compute  Va£. 

5.  If  Va£  =  0,  restart  from  step  2,  repeating  steps  3  and  4  until  the  state  and  costate 
equations  are  converged  on  the  entire  field.  Otherwise  take  =  a"  +  /(Vc,£)  and 
go  to  4. 

We  took  —  a  Va£,  where  a  is  a  parameter.  One  could  try  to  solve  the  problem 

on  the  boundary,  in  step  4  above,  using  a  direct  solver.  The  way  we  propose  here  has  the 
advantage  of  being  a  simple  modification  of  adjoint  method,  and  therefore  can  be  easily 
implemented. 

5  First  optimization  experiments 

We  introduce  a  discrete  grid  defined  as  (xi^ym)  =  {^o  +  ^A.x,7/(0)  +  mAy)  where  Ax  is 
constant  and  Ay  is  a  constant  fraction  of  the  local  height  of  the  nozzle;  see  Fig.5. 

The  steady  solution  of  the  Euler  equations  is  obtained  with  a  time- dependent  technique, 
in  the  frame  of  an  explicit  finite  volume  code.  The  conservative  variables  U  are  computed 
at  the  cell  centers,  and  the  fluxes  F  and  G  are  evaluated  at  the  cell  interfaces  using  the 
approximate  Riemann  solver  in  [14].  Higher-order  accuracy  is  achieved  using  an  Essentially 
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Non-Oscillatory  scheme  [7].  The  flow  field  values  at  cell  interfaces,  used  as  initial  conditions 
for  the  Riemann  problem,  are  reconstructed  by  means  of  a  linear  interpolation  and  using 
a  minimod  limiter.  The  amplitude  of  the  integration  step  is  chosen  according  to  the  CFL 
condition. 

The  costate  equations  are  discretized  on  the  same  grid  presented  above.  Since  they  have 
no  conservative  form,  the  numerical  solution  is  obtained  using  the  finite-difference  scheme 
proposed  in  [8]. 

The  computations  are  performed  on  a  40  x  20  grid.  Total  pressure  and  total  temperature 
at  the  inlet  are  taken  to  be  unity  and  cr(0,y)  =  0.  At  the  outlet  the  static  pressure  depends 
on  the  test-case  considered.  For  the  lower  wall  ordinate  y(0)  we  have 


y(0)  = 


0 

HUl  OCiX'‘+^  {x-\f 
0 


if  -0.5  <  a:  <  0 
if  0  <  a:  <  1 
if  1  <  a;  <  1.5 


We  try  to  recover  the  pressure  distribution  obtained  with  the  Euler  solver,  corresponding 
to  the  set  of  shape  coefficients  a  =  (0,2, 2,0).  This  means  that  the  functional  is  0  at 
the  minimum.  The  outlet  pressure  is  such  that  the  flow  presents  a  relevant  shock  at  re¬ 
compression. 

Figure  6  shows  the  functional  values  at  each  step  of  computation.  Figure  7  shows  the 
the  convergence  history  of  the  state  equation,  computed  to  second-order  accuracy,  and  the 
convergence  history  of  the  costate  equation.  Finally  in  Fig.  8,  we  present  the  starting 
pressure  distribution  and  the  one  obtained  at  the  end  of  the  optimization  procedure. 

The  practicability  of  this  approach  depends  on  the  rate  of  convergence  to  the  minimum.  In 
fact  the  state  and  costate  equations  converge  to  the  steady  solution  with  a  less  favorable  rate 
compared  to  that  of  a  simple  analysis.  It  is  easily  seen  that  since  the  shape  is  changing,  the 
flow  field  must  change  accordingly  and  so  must  the  residuals.  Figure  9  shows  a  comparison 
of  the  residuals  for  the  state  equations  in  the  case  of  a  simple  analysis  to  the  residuals  in 
the  optimization  case.  The  convergence  rate  is  3  to  4  times  slower  in  the  optimization  case. 
Considering  the  costate  equations,  the  cost  of  the  optimization  procedure  turns  to  be  of  the 
order  of  10  analyses;  using  the  first  of  the  two  algorithms  presented  in  Section  4  the  factor 
of  proportion  is  100  to  200  depending  on  the  updating  strategy  used.  The  CPU  time  needed 
on  a  DEC  3000/500  is  18  minutes  with  the  algorithm  presented.  For  the  first  algorithm  of 
Section  4,  6  hours  of  CPU  time  were  needed. 

The  present  rate  of  convergence  could  be  improved  by  changing  the  way  of  updating 
the  grid.  In  fact,  close  to  the  minimmn,  the  entire  grid  is  perturbed  to  update  only  the 
boundary.  We  believe  that,  close  to  the  minimum,  the  rate  of  convergence  can  be  improved 
by  updating  only  the  boundary  points  of  the  grid.  In  fact,  the  small  difference  between 
the  desired  pressure  and  the  one  obtained  is  due  to  the  fact  that  close  to  the  minimum 
the  convergence  rate  of  the  equations  is  reduced.  Therefore  the  pressure  p*  is  obtained 
asymptotically. 
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6  Optimal  shape  for  compressible  flows 

The  following  examples  represent  situations  for  which  the  optimal  solution  is  not  generated 
with  the  same  algorithm  we  used  to  study  the  optimization  problem.  In  the  first  case,  we 
recover  a  pressure  distribution  known  theoretically  and  compare  the  shape  obtained  with 
the  theoretical  one. 

The  Ringleb  flow  (see  [17])  is  a  two-dimensional  steady  compressible  isentropic  flow,  where 
subsonic,  transonic  and  supersonic  regions  are  represented.  It  describes  a  180°-turn  of  a 
compressible  flow;  all  the  exact  values  of  the  flow  properties  are  given  by  simple  formulas 
dependent  on  the  stream  function  and  on  the  Mach  number.  We  consider  the  portion  of 
the  flow  confined  between  two  streamlines,  which  may  be  regarded  as  solid  walls.  The 
maximum  Mach  number  on  the  bottom  streamline  is  1.6  and  the  minimum  0.8.  On  the  top 
streamline  the  maximum  Mach  number  is  0.8.  The  theoretical  Mach  number  isocontours 
for  such  a  flow  are  shown  in  Fig.  10. 

The  Ringleb  pressure  distribution  on  the  bottom  wall  is  taken  as  the  desired  distribution 
that  we  want  to  achieve. 

The  lower  wall  is  described  by  the  following  parameterization: 

y(0)  =  To  -I-  (ri  -  ro)  Y]  a,  sin(i  tt  ^ 

i=l  ~  ^0 

where  tq  is  the  distance,  measured  from  the  point  of  intersection  of  the  lines  from  the  inlet 
and  the  outlet,  to  the  first  point  on  the  lower  wall  and  ri  is  the  distance  to  the  last  point. 
The  angles  6q  and  6\  are  relative  to  the  first  and  last  point  respectively,  and  are  measured 
from  the  line  from  the  inlet. 

In  principle,  if  we  try  to  recover  the  pressure  distribution  on  the  lower  wall,  the  solution 
is  out  of  the  design  space.  We  don’t  have  any  a  priori  knowledge  of  the  values  that  the 
will  assume  and  how  close  to  the  desired  pressure  we  can  get. 

In  Fig.  11  it  is  seen  that  no  visible  difference  can  be  appreciated  between  the  theoretical 
wall  shape  and  the  optimal  shape  found.  The  points  representing  the  two  solutions  do  not 
overlap  since  they  are  computed  on  two  dilferent  grids.  In  this  case  the  functional  eq.(4) 
is  6.70  •  10“®  after  500  iterations  of  steps  2.  and  3.  of  the  second  algorithm  proposed.  See 
Fig.  12.  The  CPU  time  required  for  this  case  is  about  20  minutes. 

In  the  second  case  considered,  we  are  concerned  with  a  convergent  nozzle.  The  lower  wall 
is  represented  by  a  parameterization  similar  to  that  above.  The  inlet  Mach  number  is  2.2 
and  the  grid  is  80  x  40.  The  starting  configuration  with  pressure  contourlines  is  shown  in 
Fig.  13. 

A  relevant  shock  is  present  in  the  flow  field  and  our  objective  is  to  eliminate  it  by  requiring 
a  smooth  compression  at  the  lower  wall.  The  smooth  pressure  distribution  is  not  perfectly 
attained,  as  is  seen  in  Fig.  14.  Nevertheless  the  recompression  appears  to  be  smooth  and 
the  shock  is  eliminated  from  the  flow  field  (Fig.  15).  These  results  are  obtained  after  200 
iterations  in  about  35  minutes  of  CPU.  The  functional  is  decreased  from  3.75  •  10“^  to 
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1.35  *  10  The  computation  has  been  pursued  for  2000  iterations  and  the  functional  value 
remained  unchanged. 

Finally,  an  experiment  using  8  shape  coefficients  is  performed.  In  Fig.  16  it  is  seen  that 
the  convergence  rates  of  the  state  and  costate  equations  are  not  affected.  The  functional 
minimum  is  therefore  attained  with  the  same  number  of  iterations  as  in  the  case  of  4  shape 
coefficients. 


7  Conclusions 

The  pseudo-time  method  was  applied  to  optimization  problems  governed  by  the  Euler  equa¬ 
tions  in  two  dimensions.  The  problem  of  matching  the  pressure  distribution  to  a  desired 
one  was  considered,  both  for  subsonic  and  supersonic  flows.  The  rate  of  convergence  to 
the  minimum  for  the  cases  considered  is  3  to  4  times  slower  compared  to  that  of  the  ana¬ 
lysis  problem.  Results  were  obtained  for  Ringleb  flow  and  a  shockless  recompression  case. 
The  algorithm  could  be  implemented  with  no  substantial  changes  to  existing  adjoint  based 
codes.  Numerical  results  indicate  that  the  method  converges  at  a  rate  which  is  independent 
of  the  number  of  design  variables.  The  method  offers  a  powerful  and  inexpensive  tool  for 
the  study  of  non-intuitive  configurations  for  aerodynamic  design. 
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Figure  2;  Point  A  represents  the  desired  optimum,  point  B  the  starting  configuration.  In 
the  standard  adjoint  method,  point  A  is  reached  by  following  a  narrow  path  corresponding 
to  the  intersection  of  plane  S  and  C.  At  each  step  along  A  -  B,  the  state  and  costate 
equations  are  iterated  many  times. 


Figure  3:  In  the  new  approach  a  new  path,  A  —  B,is  taken  lying  on  the  plane  T  representing 
all  solutions  to  the  design  equation  and  the  boundary  conditions  of  state  and  costate  equa¬ 
tions.  The  computational  cost  of  working  on  this  plane  is  equivalent  to  solving  a  problem 
one  space  dimension  less  than  the  original  problem.  The  solution  of  the  state  and  costate 
equations  is  achieved  only  when  point  A  is  reached. 
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Figure  4:  Domain  of  integration  and  auxiliary  boundary 
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Logarithm  of  the  functional 


Figure  6:  Logarithm  of  the  functional  versus  the  number  of  iterations  on  the  entire  field. 


Figure  7:  Convergence  history  for  state  and  costate  equations.  Logarithm  of  the  residuals. 
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Figure  8:  Target  pressure  distribution  and  optimal  one.  Starting  pressure:  constant  distri 
bution  at  0.7  reference  value. 


Figure  9:  Convergence  history  of  state  equations  in  a  simple  analysis  compared  to  the 
convergence  of  the  state  equations  in  an  optimization  procedure. 
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Figure  10:  Ringleb  flow:  Mach  number  isocontours. 
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Figure  11:  Ringleb  flow:  starting  configuration,  theoretical  solution  and  optimal  shap  : 


Figure  14:  Convergent  nozzle:  desired  pressure  distribution  and  optimal  one. 
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Figure  16:  Convergent  nozzle.  State  and  costate  convergence  history:  4  shape  coefficients 
versus  8  shape  coefficients. 
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